Mule deer fawn recruitment dynamics in an energy disturbed landscape

Abstract Wildlife population dynamics are modulated by abiotic and biotic factors, typically climate, resource availability, density‐dependent effects, and predator–prey interactions. Understanding whether and how human‐caused disturbances shape these ecological processes is helpful for the conservation and management of wildlife and their habitats within increasingly human‐dominated landscapes. However, many jurisdictions lack either long‐term longitudinal data on wildlife populations or measures of the interplay between human‐mediated disturbance, climate, and predator density. Here, we use a 50‐year time series (1962–2012) on mule deer (Odocoileus hemionus) demographics, seasonal weather, predator density, and oil and gas development patterns from the North Dakota Badlands, USA, to investigate long‐term effects of landscape‐level disturbance on mule deer fawn fall recruitment, which has declined precipitously over the last number of decades. Mule deer fawn fall recruitment in this study represents the number of fawns per female (fawn:female ratio) that survive through the summer to October. We used this fawn recruitment index to evaluate the composite effects of interannual extreme weather conditions, energy development, and predator density. We found that density‐dependent effects and harsh seasonal weather were the main drivers of fawn fall recruitment in the North Dakota Badlands. These effects were further shaped by the interaction between harsh seasonal weather and predator density (i.e., lower fawn fall recruitment when harsh weather was combined with higher predator density). Additionally, we found that fawn fall recruitment was modulated by interactions between seasonal weather and energy development (i.e., lower fawn fall recruitment when harsh weather was combined with higher density of active oil and gas wells). Interestingly, we found that the combined effect of predator density and energy development was not interactive but rather additive. Our analysis demonstrates how energy development may modulate fluctuations in mule deer fawn fall recruitment concurrent with biotic (density‐dependency, habitat, predation, woody vegetation encroachment) and abiotic (harsh seasonal weather) drivers. Density‐dependent patterns emerge, presumably due to limited quality habitat, being the primary factor influencing fall fawn recruitment in mule deer. Secondarily, stochastic weather events periodically cause dramatic declines in recruitment. And finally, the additive effects of human disturbance and predation can induce fluctuations in fawn fall recruitment. Here we make the case for using long‐term datasets for setting long‐term wildlife management goals that decision makers and the public can understand and support.

Understanding the mechanism of wildlife response to disturbance can aid the conservation of species subjected to human pressures.
Human modification of ecosystems and fragmentation of habitats yields trophic responses across ecological communities as it forces top-down and bottom-up cascading effects on predators and prey (Karakoç et al., 2018). While ecosystems are inherently dynamic, rapid anthropogenic change has forced wildlife to occupy habitats alien to their evolutionary history resulting in changes to population dynamics across species (Guiden et al., 2019). Unfortunately, there is a lack of understanding of multiple effects of human disturbance on wildlife and their ecology. This is in some part due to a lack of longterm, systematic monitoring schemes needed to detect these trends and longitudinal studies required to identify the cause of ecosystem change in response to anthropogenic activity (Graham et al., 2005;Murphy et al., 2022).
Human disturbance to wildlife can introduce novelty into ecosystem communities and alter predator-prey interactions by influencing resource availability and modifying habitat that modulates spatial overlap, such as reducing refugia for prey species (Guiden et al., 2019;Muhly et al., 2011). How a disturbance influences predator and prey communities is dependent on the nature of the behavioral response of each group. For example, prey can benefit when predators are displaced from modified landscapes (Burr, 2014). Conversely, predators can benefit from landscape modifications that increase hunt success (Berger, 2007;Dickie et al., 2020). If both predator and prey change behaviors to avoid human activity, both species can be pushed into ever-shrinking natural habitats, where prey must alter their behavior to avoid predators, e.g., increase vigilance and other predation risk avoidance behaviors (Ciuti et al., 2012;Murphy et al., 2021). These trade-offs incur a fitness cost that can be additive to other environmental constraints. For instance, harsh seasonal weather and reduced resource availability can cause additional negative impacts to survival and recruitment of prey populations (Basille et al., 2013;Murphy et al., 2021). Therefore, it is important to disentangle the drivers of prey population dynamics within modified ecosystems, specifically the interplay between human disturbance, predation, and seasonal environmental change.
Without long-term studies on the effect of disturbance to wildlife, relatively short-term consequences of long-term effects can be misidentified that may result in human-wildlife conflict or further disturbance (Ballard et al., 2001;Graham et al., 2005;Murphy et al., 2022;Woodroffe et al., 2005). Thus, up-to-date knowledge on population productivity is important to monitor long-term population health. Recruitment can be a proxy of prey population productivity because it represents the joint contribution of fecundity and offspring survival, making it one of the most variable parameters being the primary factor influencing fall fawn recruitment in mule deer. Secondarily, stochastic weather events periodically cause dramatic declines in recruitment. And finally, the additive effects of human disturbance and predation can induce fluctuations in fawn fall recruitment. Here we make the case for using long-term datasets for setting long-term wildlife management goals that decision makers and the public can understand and support. in large herbivore populations and a key factor affecting population growth rates (Gaillard et al., 1998). Many metrics exist to examine population productivity in ungulates. For instance, adult survival is often stable whereas fecundity and recruitment have greater annual variability and have been shown to be better predictors of longterm population productivity (Gaillard et al., 1998). Recruitment is most commonly affected by climate, habitat quality, competition, and predation and may be particularly sensitive to human disturbance (Vye et al., 2017;Whittaker & Lindzey, 1999). Recruitment is also of particular interest for the conservation of hunted ungulate populations, for instance, deer species, because it can dictate harvest quotas (Ciuti et al., 2015). Consequently, studies on harvested deer populations are essential for monitoring trends and that might arise from limiting factors such as predation or disturbance (Ballard et al., 2001).
By using one of the longest term surveys of wild cervids in North America  accompanied by unique highresolution spatiotemporal data on oil and gas development, we aim to disentangle the biotic (e.g., population density and predator density) and abiotic (e.g., seasonal weather and habitat quality) drivers of mule deer (Odocoileus hemionus) fawn fall recruitment in the North Dakota Badlands. Because of the unique data on oil and gas developments, we had the opportunity to unravel whether predation by coyotes and harsh seasonal weather influenced fawn fall recruitment when interacting with long-term disturbance including a period of rapidly escalating landscape modification, i.e., rapid development of hydraulic fracturing for oil and gas (NDDMR, 2021). Mule deer fawn fall recruitment in this study represents the number of fawns per female (fawn:female.ratio) that survive to October. We used this as a proxy for fawns surviving the first 4-6 months of life where risks such as abortion, malnutrition, drought, and predation are prevalent. While fawn:female ratio in fall is a strong index of interannual fawn summer survival, it is important to state that only a portion of this population will survive the overwinter period and contribute to the adult population (Forrester & Wittmer, 2013). Our a priori hypothesis was that long-term energy disturbance would negatively affect fall fawn recruitment when combined with harsh seasonal weather and/ or predation pressure, although we did not have a priori expectations on the effect size of such interactions. Our goal was to describe the long-term effects of energy development on herbivore ecology within a landscape with increasing human-developed disturbances.

| ME THODS
We used generalized linear mixed-effect additive models (Wood, 2017) to explain the variability of mule deer fawn recruitment (fall fawn:female ratios) as a function of environmental (weather, landscape), ecological (predator, habitat), and anthropogenic disturbance (energy development, associated infrastructure) variables across 26 study sites and 50 years of monitoring in the North Dakota Badlands. We used the psych package (Revelle, 2021) in R to calculate Pearson's correlation coefficients for all predictor combinations, and we wished to use in the model and excluded collinear smooth terms (Dormann et al., 2013;Wood, 2017  Dakota's climate is continental, and temperatures are highly variable both seasonally and daily. Temperatures range from a mean daily high of 30.8°C in July (1962-2012: range 23.5-34.5°C) to a mean daily low of −16.1°C in January (range −24.9 to −6.5°C).

| Study area
Snow cover typically occurs between November and April. Mule deer share this landscape with low-density populations of pronghorn (Antilocapra americana), white-tailed deer (Odocoileus virginianus), elk (Cervus elaphus), and bighorn sheep (Ovis canadensis).
Coyotes (Canis latrans) are the primary predator of mule deer and their fawns in our study sites (Jensen, 1988). Year and mule deer density were collinear in this dataset (Pearson's correlation coefficient r p = .75, Figure 3b). Increasing deer densities over time combined with a steady decrease of fawn recruitment along the time series ( Figure 3a) suggested a densitydependent effect in the preliminary stage of our analysis. The collinearity between year and mule deer density prohibited running models with both predictors (Dormann et al., 2013). However, it was important to include the time-series (year) in the model to reduce residual temporal autocorrelation and mule deer density to assess the long-term effect of density-dependent effects on recruitment. Therefore, we ran a principal component analysis using the prcomp function within the stats package (R Development Core Team, 2021) to include information from both predictors.

|
We included both principal components which accounted for 100% of the variability of both predictors. PC1 was included to capture density-dependent effects (87.9% variability explained, Figure 4). PC1 increases in value as mule deer density increases over time, whereas PC2 encapsulated years where mule deer density decreased over time (12.1% variability explained, Figure 4).
We used the principal components of year and mule deer density to explain variation in recruitment that is not explained by densitydependent effects over time.

| Weather data
We split the year into 3 main seasons where weather is expected to drive fawn recruitment. We expect that unfavorable weather conditions recorded prior to birth (winter: November-March), during the birth season (spring: April-May) and immediately after the birth (summer: June-September) of fawns negatively affect fawn survival and, consequently, fawn recruitment recorded during the following October. Harsh winter conditions were expected to weaken female body condition with consequences on their reproductive success (Ciuti et al., 2015;Forrester & Wittmer, 2013). At the same time, snowy and cold spring weather can debilitate female body condition after fat reserves in females were depleted over winter (Forrester & Wittmer, 2013). Unfavorable spring weather extends winter, with snow sometimes persisting for longer periods and delaying greenup (Ciuti et al., 2015). Finally, hot and dry summers are expected to reduce food quality and availability thus reducing fawn survival (Jensen, 1988). We extracted annual values of average minimum temperature (°C), average precipitation (mm), and average snow depth (cm) data for winter (mean minimum temperature: −11.2°C; mean precipitation: 12.71 mm; mean snow depth: 63.76 mm) and spring (mean minimum temperature: 2.1°C; mean precipitation: 51.17 mm; mean snow depth: 5.49 mm), and average minimum temperature (°C) and average precipitation (mm) for summer (mean minimum temperature:  (Figure 2). This station has been shown to provide reliable weather data for the region previously used in mule deer population modeling (Ciuti et al., 2015).

F I G U R E 2
Map of our study area in North Dakota, showing (a) the location of mule deer study sites (in yellow) in North Dakota in relation to the Little Missouri river and the location of the Medora weather station (red triangle, middle of image, and bordering river) and (b) the total cumulative location of individual active well sites between 1962 and 2012 (black dots) and major highways (orange network) in relation to our study sites.

F I G U R E 3
Long-term population metrics for mule deer over the study period . Mule deer fawn recruitment, calculated as the ratio of fawns: female in fall (a), boxplots representing the variability among the 26 study sites (b) and mule deer density calculated as the total number of animals in spring (c) across 26 study sites in south-west North Dakota.  1962 1968 1974 1980 1986 1992 1998  To describe the heterogeneity of the environment across the 26 study sites spread along the Little Missouri River (Figure 2a), we con-  (Ciuti et al., 2014).

| Disturbance data
The density and activity of active oil and gas wells (active wells/10 km 2 ) inside mule deer survey units and within a buffer of 1 km was computed with data from satellite imagery, governmental data, and information from extraction companies for the period 1962-2012 (Figures 2b and 6). Contrary to usual energy development maps that are static over time (i.e., location of the wells drilled across a given study period), we rebuilt the life history of every well drilled from 1962 to 2012 in the mule deer study sites (Figures 2b   and 6). In our final database, we associated the density of active wells in a given year to the corresponding mule deer study site. If wells were abandoned over time (inactive), then they were not considered in the well density estimates because we assumed the level of disturbance to wildlife is not obtrusive for inactive wells. We also included road density in our models (computed as the kilometers of roads per km 2 ). This was calculated for each year in our time series  as the km of roads per km 2 of the study area to account for annual variation in road development. Annual road density in our study area was positively correlated (but not collinear) with annual active well density (Pearson's correlation coefficient r p = .19).

| Predator data
Coyote observations were collected during mule deer aerial surveys in spring as incidental observations. These coyote observations are likely the least robust data in our dataset, being derived from opportunistic observations during mule deer aerial surveys. That said, these data are unique for such a long time series, and we carefully evaluated the data prior to including it in our recruitment model.
During preliminary analysis, we considered using a relative density estimate computed as coyotes per study site; however, these did not provide robust estimates due to limited sample sizes per study site compared to the whole region. Rather, we calculated overall coyote density for the 26 study areas (hereafter also referred to as coyote density) over the course of the study period to have a consistent index of annual coyote density in spring. We used these data with the assumption that such estimates underestimate the actual coyote density but provide a long-term repeatedly collected index of interannual variation in predator presence and pressure on mule deer populations across the study area. Our coyote density estimates  were significantly positively correlated (Pearson's correlation coefficient r p = .25) with another independent index of coyote density across the Badlands derived from incidental observations by rural mail carriers (available from 1980 to 2012).

| Generalized additive mixed-effect models
We modeled the data using generalized additive mixed-effect models (Wood, 2017) using the gam function from the mgcv package in R (Wood, 2017) to examine the relationship between mule deer fawn recruitment and our selected covariates. Our final database con- are joined together at specified locations known as knots (controlled by a value K in the model). We conservatively set the K value for each covariate in the model to define the number of knots in each covariate spline, thus allowing the model to learn from the data and penalize the regression spline accordingly to avoid overfitting the data (Wood, 2006). We also included Observer and StudyArea as random effects. We opted to use a Gaussian distribution of errors (identity link) that met the assumptions of homogeneity and normality (Wood, 2017). We screened the predictor data to check for gaps in the data space coverage, which highlighted a number of variables that needed to be log10 transformed (i.e., spring precipitation, spring mule deer density, coyote density, road density, and well density). Gaps in the data space coverage can bias model estimates in generalized additive models, and thus, we followed the step-by-step protocol defined by Zuur et al. (2009), which recommends log10 transforming data to ensure good coverage. We fit alternative a priori models including one seasonal weather predictor per model (i.e., one of-average minimum temperature, average precipitation and where relevant, average snow depth-for winter, spring, and summer, respectively). Thus, we had 3 weather predictors for winter and spring and 2 for summer (where snow depth is irrelevant), resulting in a final set of 18 unique alternative models (3 winter variables x 3 spring variables x 2 summer variables). We used second-order Akaike Information Criterion (AICc) to rank alternative models and the null model. AICc values for each of the models were extracted from the Mumin package in R ( Barton, 2022); however, we did not use the Mumin package to rank the models, rather we followed classical approach to theoretic information criterion model ranking (Burnham et al., 2011).
Each model had the same a priori structure with the exception of alternative weather predictors. We designed this structure based on our understanding of the local animal ecology to examine the effects of weather, disturbances, and predators on mule deer recruitment over the 50-year period. Specifically, we modeled fawn:female ratio as a function of single terms-two principal components computed using year and mule deer density data, three principal components computed using landscape data (elevation variance, slope variance, woody vegetation percentage, longitude, and latitude), annual road density, and spring mule deer density as smoothing splines. We included interactions for seasonal weather and coyote density-to test whether predator pressure had a differential effect depending on weather harshness-and seasonal weather and well density-to Fawn: female ratio ∼ intercept + f 1 1 + … +f 13 13 + f 14 (Observer) + f 15 (StudyArea) + error (Gaussian) test whether energy development interacted with weather harshness in affecting fawn recruitment. Finally, we included an interaction between well density and coyote density to test our hypothesis of the differential ability of coyotes to prey upon fawns based on energy development-related habitat modification and loss. We also included crossed random intercepts for the observer on the plane and survey block (ID). The models were fit using the select = TRUE implementation in the gam algorithm of the mgcv package: this option adds a penalty to each smoothing term, allowing it to be penalized out of the model (when judged irrelevant by the algorithm) via optimization of the smoothing parameter selection criterion (Wood, 2017). Normality and homogeneity of GAMM residuals were successfully met by inspecting quantile-quantile (QQ) plots and simulated residual plots from the DHARMa package in R (Hartig, 2021).

F I G U R E 6
Fan chart of percentiles of oil and gas active well density (in shaded 10-percentile intervals) across 26 study areas in the North Dakota badlands for the years of record . The median density is shown as a white line. Note: Models have been ranked using the corrected Akaike Information Criterion AICc (LogLik: log-likelihood; ΔAICc: difference in AIC between any model and the top ranked one). All models had the same a priori structure (see text for full details), except for different combinations of the average seasonal weather proxies in the 3 main seasons of the mule deer annual biological cycle. We examined the model residuals for spatial (R package: sp; function: variogram) and temporal (R package: stats; function: acf) autocorrelation (Cressie, 1993;Venables & Ripley, 2002).

| RE SULTS
Of the 18 a priori generalized additive mixed-effect models fitted and compared with AICc, we selected the top-ranked model (pseudo-R 2 = 0.452; Deviance Explained = 47.8%, Table 1) containing average minimum temperature in winter, average snow depth in spring, and average minimum temperature in summer as proxies for seasonal weather. The second-ranked model was 1.18 AICc points lower than the top-ranked model (differing from the top-model only by having minimum temperature instead of snow depth as a spring weather proxy; we present the effects of spring minimum temperature in Appendix S1: Supplementary Material 2). Note that this modeling strategy has been adopted to identify the best weather indices, which were suggested by the lowest AICc value of the topranked model. We were not interested in model comparisons and averaging because structure was identical among models with the only exception being the included proxies for seasonal weather. We report approximate significance of smooth terms (including crossed random intercepts for survey observer and study area) for the topranked GAMM in Table 2. We did not find evidence of temporal (Appendix S1: Supplementary Material 3) or spatial (Appendix S1: Supplementary Material 4) autocorrelation in the residuals of the top-ranked model.
For the top-ranked model ( Table 2; see Appendix S1: Supplementary Material 5 for the smoothing splines not included as figures here), we found a negative relationship between PC1 computed using year and mule deer density, denoting a densitydependent effect (Table 2; Figure 7a). Likewise, we found PC2 for year and mule deer density to be negatively related to recruitment, this time explaining the decline in recruitment over time that was not due to density-dependent effects (Table 2; Figure 7b).
We found no effect of PC1 computed using landscape data in the model ( Table 2; Appendix S1: Supplementary Material 5A). We did find a significant effect of PC2 on recruitment (Table 2, Figure 7c); this describes an increase in recruitment in the north-eastern survey blocks less affected by woody vegetation encroachment. Finally, PC3 provided no effect on mule deer fawn recruitment ( Table 2; Appendix S1: Supplementary Material 5B). Road density also had no effect on mule deer fawn recruitment ( Table 2; Appendix S1: Supplementary Material 5C).
Active well density had a subtle but significant interaction with winter and spring seasonal weather in finely modulating mule deer recruitment. We found that active well density interacted with average minimum winter temperature (Table 2; Figure 8a), when we found higher recruitment when active wells occurred at lower densities and winters were milder (Figure 8a). Similarly, we found a significant interaction between active well density and snow depth in  Figure 8b), with snow depth strongly and inversely correlated to recruitment (particularly when exceeding ~25-30 cm of snow on the ground) while being slightly negatively impacted by active well density (Figure 8b). We found no significant effect of active well density interacting with average minimum summer temperature ( Table 2, Appendix S1: Supplementary Material 5D).
We found an interaction between coyote density and winter harshness (Table 2; Figure 9a). Fawn recruitment was higher during milder winters and when coyotes occurred at lower densities, the opposite effect on recruitment being recorded during colder winters and higher coyote densities (Figure 9a). Similarly, higher coyote densities and deeper snow depth in spring had a negative interactive effect on mule deer recruitment (Table 2, Figure 9b). We also found that hot summer temperatures, in combination with higher coyote density, led to a reduction in recruitment (Table 2, Figure 9c), with high predator density particularly impacting recruitment during hot summers.
Finally, we did not find a significant interaction between coyote density and active well density in affecting mule deer recruitment (Table 2; Figure 10). Our data suggests that coyote density and active well density are not interactive factors in this predator-prey system.
We indeed found a weak but significant negative effect of well density on fawn recruitment irrespective of coyote density. Likewise, we found a stronger negative effect of coyote density on fawn recruitment irrespective of active well density ( Figure 10). Thus, the effects of predator density and energy development appear to be additive (rather than interactive) in shaping mule deer fawn recruitment fluctuations in the North Dakota Badlands. We have included all interaction plots as single panel plots, with 95% confidence intervals shown in Appendix S1: Figure S1.

| DISCUSS ION
During the late 1800s, mule deer populations, like other big game species in the western United States, were severely reduced by over exploitation from unregulated hunting. Following the establishment of regulated hunting seasons, mule deer populations increased rapidly and peaked in the 1960s (Knue, 1991). In the early to mid-1970s, mule deer numbers began to decline again (Jensen et al., 2023). Biologistshave attributed these declines to a variety of F I G U R E 7 Effect of year and mule deer density principal components on mule deer fawn recruitment. PC1 increases (and recruitment decreases) when mule deer density increases with year of the study (a). PC2 increases (and recruitment decreases accordingly) when mule deer density decreases, accounting for years of decline not attributable to density-dependent effects (b) and the effect of landscape increasing values of PC2 on mule deer fawn recruitment, where recruitment increases in their core range (north-east) that is unaffected by the encroachment of woody vegetation such as Rocky Mountain juniper (c).  (Doyle et al., 1993;Ziegler & Myers, 1983). In contrast, black-tailed deer (O. h. columbianus and O. h. sitkensis) populations seemed less affected by these factors (Gill, 1990). Understanding the primary factors driving mule deer and other ungulate populations is essential for developing informed and effective management strategies (Forrester & Wittmer, 2013;Peek et al., 2002;Sawyer et al., 2017). oil and gas development might have a strong influence on fawn recruitment. This was not the case, rather, we found that mule deer recruitment trends were largely driven by biotic (population dynamics, habitat quality) and abiotic factors (seasonal weather severity and landscape variability) with human disturbance from oil and gas development and predation seeming to play a subtle and additive role in influencing annual oscillations in fawn recruitment.
The strong influence of density dependence on fall fawn recruitment over time suggests limited quality habitat to be the primary driver of mule deer in this northern Great Plains population, with Rocky Mountain juniper encroachment being a key factor of habitat degradation for mule deer. Density dependence is further supported by the observation that within Theodore Roosevelt National Park, a nonhunted population with a mule deer density 69% higher than outside the park, had a fall fawn recruitment rate that was 24% lower than through the Badlands (Jensen et al., 2022). Density-dependent effects in wildlife are manifested in a series of feedback mechanisms that lead to fluctuations in populations depending on the relative density of the animal. Density-dependent effects are well studied in wildlife (Hixon & Webster, 2002;Tanner, 1966) and, in particular, in wild ungulate populations (Brown, 2011;Gaillard et al., 1998;Saether, 1997). Here we found another example of densitydependent effects in North Dakota mule deer, where the effect likely modulates populations through intraspecific competition for high-quality habitats. We also observed greater declines in recruitment in areas where primary mule deer habitat has been affected by F I G U R E 9 Interactive effects of coyote density and average winter minimum temperature (a), coyote density and average snow depth in spring (b), and coyote density and average summer minimum temperature (c) on fall mule deer fawn recruitment (fawn/female ratio) in the North Dakota badlands  as predicted by the top ranked generalized additive mixed effect model. condition is closely related to the quality of available resources, in contrast to high quantities of less nutritious foods (Wallmo, 1981).
Mule deer have been observed to adjust their dietary niche when populations occur at higher densities with evidence of niche expansion (Stewart et al., 2011) and contraction (Nicholson et al., 2006), resulting in negative consequences for population productivity.
Deer use fat and protein reserves accumulated during late summer and fall to survive the winter. As a result, behavioral changes due to high population density, or trade-offs due to perceived predation risk, have fitness and reproductive consequences due to suboptimal body condition. Effects of reduced body condition are exacerbated by harsh seasonal weather, habitat integrity, and intraspecific competition for high-quality forage (Christie et al., 2015;Ciuti et al., 2015;Gamo & Beck, 2017).
Following density dependence, stochastic severe winter weather, particularly cold temperatures, and snow depths >25 cm were found to have secondary impacts on recruitment of fawns to the following fall. Harsh winter weather is known to weaken female body condition and decreases reproductive success, and unfavorable spring weather prolongs negative impacts of winter conditions, for instance, when snow persists and delays green-up of foraging sites (Ciuti et al., 2015;Parker et al., 1996;. Whereas, hot and dry summers were expected to reduce food quality and availability, resulting in reduced fawn survival and subsequently lower autumn recruitment (Ciuti et al., 2015;Hurley et al., 2011). Seasonal weather affects energy expenditure (Freddy et al., 1986) and forage quality (Laycoek & Price, 1970) and subsequently body condition . Hot summers and higher coyote densities were found to accompany lower fawn recruitment. Drought conditions could have resulted in poor body condition and less vegetative concealment cover for fawns and thus increased fawn susceptibility to predation (Ballard et al., 2001). Gullikson (2019) reported avoidance of oil wells by white-tailed deer (Odocoileus virginianus) in western North Dakota, but physical barriers or suitable concealment cover habitat adjacent to well pads can allow deer to tolerate wells and associated activity.
Despite predators, coyotes in our study, being a primary source of direct mortality for both adults and juveniles (Forrester & Wittmer, 2013;Lingle et al., 2008), it is rare that predation causes long-term population declines in mule deer (Forrester & Wittmer, 2013). A review of deer-predator relationships (Ballard et al., 2001) reported that coyotes, mountain lions (Puma concolor), and wolves might be substantial mortality factors in some areas under certain conditions, but large herbivores at, or near, carrying capacity do not exhibit strong responses to predator removal, as emphasized recently by Bowyer et al. (2005Bowyer et al. ( , 2014Bowyer et al. ( , 2020. In our study, in years where we recorded harsh seasonal weather conditions, coyotes had a measurable effect on reducing fall mule deer recruitment. Coyote-induced mortality to mule deer populations is typically compensatory when prey density is high (i.e., when predator-induced mortality increases, other forms of mortality typically decrease, Bergman et al., 2015), therefore it will not reduce populations. Coyote predation can be additive (i.e., capable of limiting recruitment supplementary to other forms of mortality) if other sources of proximate mortality are higher, such as weather severity (Bergman et al., 2015). In addition to direct predation of mule deer fawns, coyotes may also influence mule deer behavior and space use through nonconsumptive effects (Laundré et al., 2010), also called a trait-mediated effect which cause behavioral changes, often with fitness consequences, as prey balance the need for high-quality forage, protective habitat, and safety from perceived predation risk (Frid & Dill, 2002;Preisser & Bolnick, 2008;Walther, 1969). We acknowledge that the data collected on coyotes for this study were incidental observations during mule deer surveys: that said, these data provide a good proxy for overall coyote trends over the course F I G U R E 1 0 Interactive effect of coyote density and well density on fall mule deer fawn recruitment (fawn/female ratio) in the North Dakota badlands  as predicted by the top ranked generalized additive mixed effect model. Note that this interaction was not significant in the top-ranked model. of this 50-year study; however, they likely do not accurately quantify the exact impact of predation on mule deer recruitment across a larger region. We therefore recommended that a targeted data collection of coyote population trends across mule deer jurisdictional ranges needed to quantify the degree to which coyote predation causes additive mortality for mule deer fawns.
Our data suggest that human disturbance might play a similar role in modulating mule deer fitness if the placement of oil and gas wells alters mule deer habitats and behavior. We observed, for the first time using high-quality spatial data depicting the life-history of all wells in the region over 50 years, that there was an effect of active well density which was significant when interacting with average minimum winter temperature and with average snow depth in spring. Human disturbance and the modification of ecosystems with human infrastructure can change prey behavior, as it is possible for prey to respond to human activity in the same way they respond to predators (Ellrich et al., 2015;Frid & Dill, 2002;Kolar et al., 2017;Sawyer et al., 2017). This phenomenon has been observed in many marine and terrestrial species displaying avoidance of human presence and related infrastructures (Allen & Read, 2000;de la Torre et al., 2000;Dyer et al., 2001) as well as in mule deer avoiding energy development sites across their range (Kolar et al., 2017;Lynch et al., 2015;Sawyer et al., 2017). Thus, in years where winter weather is harsh and persists into spring, the effect of landscape modification from energy development appears to impact mule deer fitness but does not alter predator-prey relationships. Contrary to our expectations, and as suggested by earlier analysis (Ciuti et al., 2014), our models did not detect a relationship between high densities of active oil and gas wells and alterations to the success of coyote predation on mule deer fawns, even at high predator density. Although we did not detect an effect of roads directly on recruitment in our models, it is likely that human modification of the landscape resulting from energy development shapes predator-prey relationships by influencing habitat quality and availability and through the introduction of novel landscape features such as roads and well pads (DeMars & Boutin, 2018;Guiden et al., 2019). In essence, given a year where proximate causes of mortality are high (e.g., high population density, low available habitat quality, and harsh seasonal weather) and the population is under duress, the contribution of high predator density and/or human disturbance appears to play a role in further limiting populations in that year but this does not drive the observed longterm decline in fawn fall recruitment.
Documenting the relative importance of factors on population trends in ungulates can be a complex task with interactions among covariates (e.g., population density, habitat availability, climate, predators, and human disturbance) being difficult to disentangle (Bergman et al., 2015;Forrester & Wittmer, 2013;Gaillard et al., 2000;Peek, 1980;Sinclair & Krebs, 2002). Here we demonstrate the value of long-term systematic wildlife monitoring datasets for understanding multivariate effects on wildlife population dynamics in a changing landscape. Quantifying and understanding the population-level repercussions of disturbances and other contributing factors on animal life history, behavior, and physiology is a conservation priority across the world not unique to energy development. It is a requirement in many fields, such as game management, sustainable development, and as part of Environmental Impact Assessment (EIA), which is required by law in the European Union (European Habitats Directive 92/43/EEC) and United States (National Environmental Protection Act, 42 U.S.C. § § 4321 et seq; Marine Mammal Protection Act, 16 U.S.C. § § 1361 et seq.) legislation (King et al., 2015;Pirotta et al., 2018). Unfortunately, comprehensive risk assessments that encompass multivariate stressors on wildlife populations over time are rarely available due to a lack of long-term data collection and lack of partnerships that facilitate the appropriate analysis of such data. This dataset identifies the long-term relative importance and ranking of various factors that primarily drive this mule deer fawn population that is availability of quality habitat and severe weather, with human disturbance and predation playing lesser roles. Chronic wasting disease may soon become another primary driver of mule deer populations (Conner et al., 2021;Gillin & Mawdsley, 2018;Potapov et al., 2016).
Our paper utilizes mule deer fawn to female ratios (fawn:female) in fall as an index of fawns surviving the hazardous first 4-6 months of life where mortality causes such as abortion, malnutrition, drought, and predation must be confronted and survived. Attempts at predator control have largely been shown to be ineffective, costly and are increasingly being negatively viewed by the public (Ballard et al., 2001;Jensen et al., 2023). An important consideration here is the correlation between mule deer fawn fall recruitment (prewinter) and post-winter survival relative to how it contributes to the adult population. Future research should focus on correlating overwinter fawn survival to fall fawn recruitment to allow wildlife managers to robustly use this index to predict future adult population recruitment. Rapidly changing ecological conditions and landscapes will inevitably make wildlife management even more difficult it would seem, therefore, that we need to follow a data-driven approach to wildlife conservation and management, allowing robust analysis of long-term monitoring data to explain the science in an empirical manner, and broaden our base of support for sound management (Prot, 2015). Management agencies and decision makers need to focus on aspects that can be controlled and have long-term effects, such as maintaining quality habitat and mitigating human disturbance. As human populations continue to grow and develop supporting infrastructure (e.g., roads, housing and energy developments, and conversion of land from native vegetation to cropland) within and around wildlife habitat, we recognize the value of population monitoring schemes and the use of methodologies to understand the consequences of human development.

ACK N OWLED G M ENTS
Long-term datasets of big game population demographics are rare,